4-EM

Plot of the data

library('mvtnorm')

plot.mixture <- function(model, YY, ...) {
  plot(YY, ...)
  mmu <- model$mmu
  SS <- model$SS
  m <- length(mmu)
  xx <- seq(min(YY[,1]), max(YY[,1]), by=.1)
  yy <- seq(min(YY[,2]), max(YY[,2]), by=.1)
  for (ell in 1:m) {
    contour(xx, yy, 
             outer(xx, yy, function(xx, yy) dmvnorm(cbind(xx, yy), mmu[[ell]], SS[[ell]])),
             add=T, col=1+ell%%7)
  }
}

YY <- as.matrix(faithful)
par(mfrow=c(1,2))
# boxplot(YY$eruptions)
# boxplot(YY$waiting)
par(mfrow=c(1,1))
plot(YY)

EM & LLK

llk <- function(model, YY) {
  mmu <- model$mmu
  SS <- model$SS
  pp <- model$pp
  m <- length(mmu)
  N <- nrow(YY)
  PPHI <- matrix(0, N, m)
  for (ell in 1:m) {
    PPHI[,ell] <- dmvnorm(YY, mean=mmu[[ell]], sigma=SS[[ell]])
  }
  sum(log(PPHI %*% pp))
}

em.gauss <- function(m, YY, llk.diff.eps=0.1, do.plot=FALSE, verbose=FALSE) {
  d <- ncol(YY)
  N <- nrow(YY)
  phi <- dmvnorm
  mmu <- list()
  SS <- list()
  ww <- list()
  zz <- list()
  llk.prev <- -Inf
  llk.hist <- NULL
  # Initialization
  # TODO: Split the data into m clusters at random; calculate mmu, SS, pp for each cluster
  for (ell in 1:m) {
    # pick two data points at random
    mmu[[ell]] <- sample(as.data.frame(t(YY)), 1)[[1]]
    SS[[ell]] <- cov(YY)
  }
  pp <- rep(1/m, m)
  if (do.plot) {
    plot.mixture(list(pp=pp, mmu=mmu, SS=SS), YY, main=paste('n.iter =', 0))
  }
  n.iter <- 1
  repeat {
    # Expectation
    PPHI <- matrix(0, N, m)
    for (ell in 1:m) {
      PPHI[,ell] <- phi(YY, mean=mmu[[ell]], sigma=SS[[ell]])
    }
    for (ell in 1:m) {
      zz[[ell]] <- pp[[ell]] * PPHI[,ell]  / PPHI %*% pp
    }
    for (ell in 1:m) {
      ww[[ell]] <- zz[[ell]] / sum(zz[[ell]])
    }
    # Maximization
    for (ell in 1:m) {
      mmu[[ell]] <- as.vector(t(YY) %*% ww[[ell]])
      pp[[ell]] <- sum(zz[[ell]]) / N
      # FIXME: vectorize
      for (i1 in 1:d) {
        for (i2 in 1:d) {
          SS[[ell]][i1, i2] <- sum(ww[[ell]] * (YY[, i1] - mmu[[ell]][i1]) * (YY[, i2] - mmu[[ell]][i2]))
        }
      }
    }
    llk.cur <- llk(list(pp=pp, mmu=mmu, SS=SS), YY)
    llk.diff <- abs(llk.cur - llk.prev)
    if (verbose) {
      cat('# ', n.iter, ' LLK = ', llk.cur, ' (', llk.diff, ')\n', sep='')
    }
    llk.hist <- c(llk.hist, llk.cur)
    if (do.plot) {
      plot.mixture(list(pp=pp, mmu=mmu, SS=SS), YY, main=paste('n.iter =', n.iter))
    }
    if (llk.diff < llk.diff.eps) {
      if (do.plot) {
        plot(llk.hist, type='o', col='blue', main='Log likelihood', xlab='n.iter')
      }
      return(list(pp=pp, mmu=mmu, SS=SS, zz=zz, llk=llk.cur, llk.hist=llk.hist))
    } else {
      llk.prev <- llk.cur
      n.iter <- n.iter + 1
    }
  }
}

EM for \(m = 2\)

set.seed(2)
model.2 <- em.gauss(2, YY, do.plot=TRUE, verbose=TRUE)

## # 1 LLK = -1266.242 (Inf)

## # 2 LLK = -1235.165 (31.0769)

## # 3 LLK = -1207.979 (27.1865)

## # 4 LLK = -1200.716 (7.263274)

## # 5 LLK = -1196.433 (4.282454)

## # 6 LLK = -1192.152 (4.281588)

## # 7 LLK = -1186.858 (5.294085)

## # 8 LLK = -1179.499 (7.358276)

## # 9 LLK = -1168.903 (10.59576)

## # 10 LLK = -1154.545 (14.3581)

## # 11 LLK = -1140.924 (13.62177)

## # 12 LLK = -1131.893 (9.030789)

## # 13 LLK = -1130.32 (1.572512)

## # 14 LLK = -1130.267 (0.05378514)

model.2$pp

## [1] 0.6438597 0.3561403

model.2$mmu

## [[1]]
## [1]  4.290234 79.974944
## 
## [[2]]
## [1]  2.037046 54.485312

model.2$SS

## [[1]]
##           eruptions    waiting
## eruptions 0.1692466  0.9315412
## waiting   0.9315412 35.9457512
## 
## [[2]]
##            eruptions   waiting
## eruptions 0.06969635  0.440863
## waiting   0.44086299 33.739059

EM for \(m = 3\)

set.seed(1)
model.3 <- em.gauss(3, YY, do.plot=TRUE, verbose=TRUE)

## # 1 LLK = -1286.069 (Inf)

## # 2 LLK = -1285.762 (0.3072853)

## # 3 LLK = -1285.387 (0.3742914)

## # 4 LLK = -1284.823 (0.5646781)

## # 5 LLK = -1283.835 (0.9871156)

## # 6 LLK = -1281.761 (2.074725)

## # 7 LLK = -1276.226 (5.534951)

## # 8 LLK = -1258.127 (18.09884)

## # 9 LLK = -1219.377 (38.75022)

## # 10 LLK = -1199.715 (19.66186)

## # 11 LLK = -1188.31 (11.40481)

## # 12 LLK = -1175.743 (12.56663)

## # 13 LLK = -1161.195 (14.54839)

## # 14 LLK = -1145.295 (15.8997)

## # 15 LLK = -1132.672 (12.62365)

## # 16 LLK = -1127.521 (5.15066)

## # 17 LLK = -1123.44 (4.080583)

## # 18 LLK = -1121.309 (2.131676)

## # 19 LLK = -1120.531 (0.7774306)

## # 20 LLK = -1120.231 (0.3000882)

## # 21 LLK = -1120.087 (0.1446344)

## # 22 LLK = -1119.999 (0.08806301)

Fit for several initial parameters

Choose the best model out of several runs with the highest LLK.

fit.em.gauss <- function(m, YY, n.tries=10, verbose=FALSE) {
  prev.llk <- -Inf
  model.best <- NULL
  for (n.iter in 1:n.tries) {
    if (verbose) {
      cat('= Fitting model', n.iter, '...\n')
    }
    model <- em.gauss(m, YY, verbose=verbose)
    if (model$llk > prev.llk) {
      model.best <- model
    }
  }
  return(model)
}

set.seed(1)
model <- fit.em.gauss(3, YY, verbose=TRUE)

## = Fitting model 1 ...
## # 1 LLK = -1286.069 (Inf)
## # 2 LLK = -1285.762 (0.3072853)
## # 3 LLK = -1285.387 (0.3742914)
## # 4 LLK = -1284.823 (0.5646781)
## # 5 LLK = -1283.835 (0.9871156)
## # 6 LLK = -1281.761 (2.074725)
## # 7 LLK = -1276.226 (5.534951)
## # 8 LLK = -1258.127 (18.09884)
## # 9 LLK = -1219.377 (38.75022)
## # 10 LLK = -1199.715 (19.66186)
## # 11 LLK = -1188.31 (11.40481)
## # 12 LLK = -1175.743 (12.56663)
## # 13 LLK = -1161.195 (14.54839)
## # 14 LLK = -1145.295 (15.8997)
## # 15 LLK = -1132.672 (12.62365)
## # 16 LLK = -1127.521 (5.15066)
## # 17 LLK = -1123.44 (4.080583)
## # 18 LLK = -1121.309 (2.131676)
## # 19 LLK = -1120.531 (0.7774306)
## # 20 LLK = -1120.231 (0.3000882)
## # 21 LLK = -1120.087 (0.1446344)
## # 22 LLK = -1119.999 (0.08806301)
## = Fitting model 2 ...
## # 1 LLK = -1234.707 (Inf)
## # 2 LLK = -1154.135 (80.57213)
## # 3 LLK = -1128.975 (25.15987)
## # 4 LLK = -1127.784 (1.190921)
## # 5 LLK = -1126.533 (1.251448)
## # 6 LLK = -1125.433 (1.099579)
## # 7 LLK = -1124.573 (0.8599739)
## # 8 LLK = -1123.894 (0.678957)
## # 9 LLK = -1123.329 (0.565005)
## # 10 LLK = -1122.858 (0.4711058)
## # 11 LLK = -1122.491 (0.3669589)
## # 12 LLK = -1122.223 (0.2675441)
## # 13 LLK = -1122.026 (0.1975011)
## # 14 LLK = -1121.868 (0.1579442)
## # 15 LLK = -1121.73 (0.137739)
## # 16 LLK = -1121.603 (0.1276629)
## # 17 LLK = -1121.48 (0.1225137)
## # 18 LLK = -1121.36 (0.1196227)
## # 19 LLK = -1121.243 (0.1176387)
## # 20 LLK = -1121.127 (0.1158781)
## # 21 LLK = -1121.013 (0.1140001)
## # 22 LLK = -1120.901 (0.111843)
## # 23 LLK = -1120.792 (0.1093402)
## # 24 LLK = -1120.685 (0.1064765)
## # 25 LLK = -1120.582 (0.1032649)
## # 26 LLK = -1120.482 (0.09973455)
## = Fitting model 3 ...
## # 1 LLK = -1277.145 (Inf)
## # 2 LLK = -1259.889 (17.25532)
## # 3 LLK = -1218.301 (41.58858)
## # 4 LLK = -1179.961 (38.33938)
## # 5 LLK = -1145.735 (34.2264)
## # 6 LLK = -1124.045 (21.68996)
## # 7 LLK = -1121.814 (2.230887)
## # 8 LLK = -1120.946 (0.8686486)
## # 9 LLK = -1120.557 (0.3884166)
## # 10 LLK = -1120.364 (0.1927393)
## # 11 LLK = -1120.24 (0.1243156)
## # 12 LLK = -1120.142 (0.09763961)
## = Fitting model 4 ...
## # 1 LLK = -1272.919 (Inf)
## # 2 LLK = -1242.463 (30.45639)
## # 3 LLK = -1202.45 (40.01349)
## # 4 LLK = -1183.286 (19.16338)
## # 5 LLK = -1161.098 (22.1881)
## # 6 LLK = -1133.807 (27.29108)
## # 7 LLK = -1125.533 (8.27415)
## # 8 LLK = -1122.386 (3.146323)
## # 9 LLK = -1120.903 (1.483316)
## # 10 LLK = -1120.117 (0.7864931)
## # 11 LLK = -1119.695 (0.42118)
## # 12 LLK = -1119.475 (0.2209002)
## # 13 LLK = -1119.36 (0.1149018)
## # 14 LLK = -1119.299 (0.06034246)
## = Fitting model 5 ...
## # 1 LLK = -1288.577 (Inf)
## # 2 LLK = -1288.068 (0.5084829)
## # 3 LLK = -1287.597 (0.4707323)
## # 4 LLK = -1286.995 (0.6021216)
## # 5 LLK = -1286.026 (0.9691306)
## # 6 LLK = -1284.114 (1.911752)
## # 7 LLK = -1279.376 (4.738229)
## # 8 LLK = -1264.716 (14.66069)
## # 9 LLK = -1228.145 (36.57005)
## # 10 LLK = -1206.73 (21.41557)
## # 11 LLK = -1199.818 (6.911756)
## # 12 LLK = -1194.515 (5.303474)
## # 13 LLK = -1187.251 (7.263575)
## # 14 LLK = -1174.669 (12.58193)
## # 15 LLK = -1153.767 (20.90242)
## # 16 LLK = -1131.888 (21.87843)
## # 17 LLK = -1123.439 (8.449541)
## # 18 LLK = -1120.886 (2.552951)
## # 19 LLK = -1120.271 (0.6144042)
## # 20 LLK = -1120.081 (0.1908042)
## # 21 LLK = -1119.989 (0.09145503)
## = Fitting model 6 ...
## # 1 LLK = -1288.16 (Inf)
## # 2 LLK = -1287.846 (0.3134081)
## # 3 LLK = -1287.449 (0.397249)
## # 4 LLK = -1286.875 (0.5743079)
## # 5 LLK = -1285.943 (0.9320087)
## # 6 LLK = -1284.168 (1.774968)
## # 7 LLK = -1279.969 (4.198858)
## # 8 LLK = -1268.352 (11.61639)
## # 9 LLK = -1248.586 (19.76627)
## # 10 LLK = -1234.311 (14.27505)
## # 11 LLK = -1220.973 (13.33768)
## # 12 LLK = -1209.253 (11.72008)
## # 13 LLK = -1199.989 (9.264551)
## # 14 LLK = -1191.219 (8.769459)
## # 15 LLK = -1178.88 (12.33982)
## # 16 LLK = -1157.449 (21.43039)
## # 17 LLK = -1130.774 (26.67494)
## # 18 LLK = -1121.587 (9.186763)
## # 19 LLK = -1120.036 (1.551722)
## # 20 LLK = -1119.734 (0.3017132)
## # 21 LLK = -1119.621 (0.1131207)
## # 22 LLK = -1119.55 (0.0708383)
## = Fitting model 7 ...
## # 1 LLK = -1287.116 (Inf)
## # 2 LLK = -1286.855 (0.2609048)
## # 3 LLK = -1286.738 (0.1176217)
## # 4 LLK = -1286.673 (0.06533044)
## = Fitting model 8 ...
## # 1 LLK = -1277.105 (Inf)
## # 2 LLK = -1264.75 (12.3548)
## # 3 LLK = -1231.585 (33.16466)
## # 4 LLK = -1180.806 (50.77975)
## # 5 LLK = -1159.923 (20.88298)
## # 6 LLK = -1146.501 (13.42163)
## # 7 LLK = -1132.846 (13.65529)
## # 8 LLK = -1128.071 (4.774265)
## # 9 LLK = -1127.709 (0.3620281)
## # 10 LLK = -1127.572 (0.1376114)
## # 11 LLK = -1127.436 (0.1360656)
## # 12 LLK = -1127.292 (0.143773)
## # 13 LLK = -1127.14 (0.1523581)
## # 14 LLK = -1126.979 (0.1608453)
## # 15 LLK = -1126.81 (0.1688024)
## # 16 LLK = -1126.634 (0.1759478)
## # 17 LLK = -1126.452 (0.1821249)
## # 18 LLK = -1126.265 (0.1872763)
## # 19 LLK = -1126.073 (0.1913882)
## # 20 LLK = -1125.879 (0.1944187)
## # 21 LLK = -1125.683 (0.1962423)
## # 22 LLK = -1125.486 (0.1966388)
## # 23 LLK = -1125.291 (0.1953392)
## # 24 LLK = -1125.098 (0.1921203)
## # 25 LLK = -1124.912 (0.1869126)
## # 26 LLK = -1124.732 (0.179869)
## # 27 LLK = -1124.56 (0.1713491)
## # 28 LLK = -1124.398 (0.1618245)
## # 29 LLK = -1124.247 (0.1517568)
## # 30 LLK = -1124.105 (0.1415104)
## # 31 LLK = -1123.974 (0.1313287)
## # 32 LLK = -1123.853 (0.1213551)
## # 33 LLK = -1123.741 (0.111671)
## # 34 LLK = -1123.639 (0.1023272)
## # 35 LLK = -1123.545 (0.09336253)
## = Fitting model 9 ...
## # 1 LLK = -1233.099 (Inf)
## # 2 LLK = -1186.704 (46.39474)
## # 3 LLK = -1164.808 (21.89594)
## # 4 LLK = -1145.595 (19.21353)
## # 5 LLK = -1134.166 (11.42824)
## # 6 LLK = -1126.752 (7.414515)
## # 7 LLK = -1123.039 (3.712398)
## # 8 LLK = -1121.495 (1.544885)
## # 9 LLK = -1120.837 (0.6574851)
## # 10 LLK = -1120.508 (0.3292761)
## # 11 LLK = -1120.311 (0.1967074)
## # 12 LLK = -1120.178 (0.1327207)
## # 13 LLK = -1120.082 (0.09668214)
## = Fitting model 10 ...
## # 1 LLK = -1246.34 (Inf)
## # 2 LLK = -1201.032 (45.3079)
## # 3 LLK = -1179.25 (21.78218)
## # 4 LLK = -1161.035 (18.21422)
## # 5 LLK = -1139.702 (21.33387)
## # 6 LLK = -1126.924 (12.7777)
## # 7 LLK = -1123.16 (3.763832)
## # 8 LLK = -1121.775 (1.385447)
## # 9 LLK = -1121.157 (0.617467)
## # 10 LLK = -1120.82 (0.3373408)
## # 11 LLK = -1120.592 (0.2281166)
## # 12 LLK = -1120.413 (0.178957)
## # 13 LLK = -1120.262 (0.1510623)
## # 14 LLK = -1120.13 (0.1311716)
## # 15 LLK = -1120.016 (0.11493)
## # 16 LLK = -1119.915 (0.1009953)
## # 17 LLK = -1119.826 (0.08892939)

Model Selection

df <- function(model) {
  m <- length(model$mmu)
  d <- length(model$mmu[[1]])
  (m - 1) + m * d + m * d*(d+1) / 2
}

aic.gm <- function(model, YY) {
  2 * (llk(model, YY) - df(model))
}

bic.gm <- function(model, YY) {
  N <- nrow(YY)
  2 * llk(model, YY) - df(model) * log(N)
}

Test df, AIC, BIC

model <- em.gauss(3, YY, do.plot=FALSE)
df(model)

## [1] 17

aic.gm(model, YY)

## [1] -2273.804

bic.gm(model, YY)

## [1] -2335.103

Select the model (\(m\)) according to IC

ic.mod.sel <- function(m.max, n.tries=10, do.plot=FALSE) {
  aic <- numeric(m.max)
  bic <- numeric(m.max)
  for (m in 1:m.max) {
    cat('m =', m, '\n')
    model <- fit.em.gauss(m, YY, n.tries=n.tries, verbose=FALSE)
    # model <- em.gauss(m, YY)
    aic[m] <- aic.gm(model, YY)
    bic[m] <- bic.gm(model, YY)
  }
  if (do.plot) {
    matplot(cbind(aic, bic), type='b', pch=c(1,2), col=c('blue', 'red'), lty=c('solid'),
            main='', xlab='m')
  }
  legend("topleft", bty = "n",
         legend = c('AIC', 'BIC'),
         col = c('blue', 'red'),
         lty = c('solid', 'solid'))
  return(list(aic=aic, bic=bic))
}

set.seed(1)
ic.mod.sel(7, 20, do.plot=TRUE)

## m = 1 
## m = 2 
## m = 3 
## m = 4 
## m = 5 
## m = 6 
## m = 7

## $aic
## [1] -2589.593 -2282.528 -2273.093 -2276.686 -2265.537 -2264.694 -2276.381
## 
## $bic
## [1] -2607.623 -2322.192 -2334.392 -2359.619 -2370.105 -2390.897 -2424.218

Cluster data points according to \(\mathbf{z}\)

plot.classification <- function(model, YY, plot.densities=TRUE) {
  N <- nrow(YY)
  m <- length(model$pp)
  ZZ <- matrix(unlist(model$zz), N, m)
  xlim <- c(min(YY[,1]), max(YY[,1]))
  ylim <- c(min(YY[,2]), max(YY[,2]))
  if (plot.densities) {
    plot.mixture(model, YY, main=paste('m =', m))
  } else {
    plot(YY, type='n', main=paste('m =', m))
  }
  class.labels <- apply(ZZ, 1, which.max)
  for (ell in 1:m) {
    par(new=TRUE)
    plot(YY[class.labels==ell,], col=(1+ell), xlab='', ylab='', xlim=xlim, ylim=ylim, axes=FALSE)
  }
}

set.seed(1)
model <- fit.em.gauss(2, YY, 20, verbose=FALSE)
ZZ <- matrix(unlist(model$zz), 272, 2)
head(ZZ)

##              [,1]         [,2]
## [1,] 1.000000e+00 3.390575e-09
## [2,] 1.609446e-09 1.000000e+00
## [3,] 9.999897e-01 1.029217e-05
## [4,] 9.543642e-06 9.999905e-01
## [5,] 1.000000e+00 1.849258e-21
## [6,] 6.361447e-03 9.936386e-01

head(rowSums(ZZ))

## [1] 1 1 1 1 1 1

# matplot(ZZ)
plot.classification(model, YY)

par(mfrow=c(2,2))
for (ell in 2:5) {
  set.seed(1)
  plot.classification(fit.em.gauss(ell, YY, 30, verbose=FALSE), YY, plot.densities=FALSE)
}

par(mfcol=c(1,1))

Compare with Mclust

library('mclust')

## Package 'mclust' version 5.2.3

## Type 'citation("mclust")' for citing this R package in publications.

Classification for \(m=2\)

set.seed(1)
par(mfrow=c(2,2))
for (ell in 2:5) {
  set.seed(1)
  M <- Mclust(YY, G=ell, modelNames='VVV')
  plot(M, what='classification')
}

par(mfcol=c(1,1))

par(mfrow=c(2,2))
for (ell in 2:5) {
  set.seed(1)
  M <- Mclust(YY, G=ell, modelNames='VVV')
  plot(M, what='uncertainty')
}

par(mfcol=c(1,1))

par(mfrow=c(2,2))
for (ell in 2:5) {
  set.seed(1)
  M <- Mclust(YY, G=ell, modelNames='VVV')
  plot(M, what='density', type = 'image', col = 'dodgerblue3', grid = 100)
  # plot(M.2.VVV, what='density', type = 'persp')
}

par(mfcol=c(1,1))

df, LLK, BIC for \(m=2\)

set.seed(1)
model.2 <- fit.em.gauss(2, YY, n.tries=20, verbose=FALSE)
plot.mixture(model.2, YY)

df(model.2)

## [1] 11

# aic.gm(model, YY)
bic.gm(model.2, YY)

## [1] -2322.192

llk(model.2, YY)

## [1] -1130.264

M.2.VVV <- Mclust(YY, G=2, modelNames='VVV')
M.2.VVV$df

## [1] 11

M.2.VVV$bic

## [1] -2322.192

M.2.VVV$loglik

## [1] -1130.264

BIC for several \(m\)s

set.seed(1)
my.bic <- ic.mod.sel(7, 20, do.plot=TRUE)$bic

## m = 1 
## m = 2 
## m = 3 
## m = 4 
## m = 5 
## m = 6 
## m = 7

set.seed(1)
mc.bic <- mclustBIC(YY, modelNames='VVV')
cat('Our BICs:\n')

## Our BICs:

print(my.bic)

## [1] -2607.623 -2322.192 -2334.392 -2359.619 -2370.105 -2390.897 -2424.218

cat('Mclust\'s BICs:\n')

## Mclust's BICs:

print(as.numeric(mc.bic))

## [1] -2607.623 -2322.192 -2333.894 -2359.216 -2385.288 -2398.974 -2426.488
## [8] -2434.990 -2447.286

matplot(cbind(mc.bic[1:7], my.bic), pch=c(1,2), type='b')

\(\boldsymbol{\Sigma}\) structure according to Mclust

BIC <- mclustBIC(YY)
plot(BIC)

summary(BIC)

## Best BIC values:
##              EEE,3        EEE,4        VVE,2
## BIC      -2314.386 -2320.207047 -2320.432965
## BIC diff     0.000    -5.820699    -6.046616

BIC

## Bayesian Information Criterion (BIC):
##         EII       VII       EEI       VEI       EVI       VVI       EEE
## 1 -4024.721 -4024.721 -3055.835 -3055.835 -3055.835 -3055.835 -2607.623
## 2 -3452.998 -3458.300 -2354.601 -2350.607 -2352.618 -2346.065 -2325.220
## 3 -3377.712 -3336.542 -2323.008 -2332.698 -2332.204 -2342.371 -2314.386
## 4 -3230.246 -3245.732 -2323.676 -2331.829 -2334.756 -2343.068 -2320.207
## 5 -3149.395 -3128.214 -2337.696 -2348.300 -2355.891 -2374.307 -2336.975
## 6 -3081.411 -3067.559 -2338.118 -2363.112 -2357.725 -2372.748 -2347.297
## 7 -2990.335 -2998.497 -2356.461 -2370.061 -2375.851 -2393.101 -2361.206
## 8 -2978.090 -2991.851 -2371.809        NA -2395.989        NA -2376.917
## 9 -2899.779 -2920.951 -2388.629        NA -2399.083        NA -2393.728
##         EVE       VEE       VVE       EEV       VEV       EVV       VVV
## 1 -2607.623 -2607.623 -2607.623 -2607.623 -2607.623 -2607.623 -2607.623
## 2 -2324.273 -2322.972 -2320.433 -2329.116 -2325.416 -2327.598 -2322.192
## 3 -2322.690 -2322.094 -2332.433 -2338.986 -2329.352 -2335.618 -2333.894
## 4 -2340.044 -2332.885 -2354.890 -2336.750 -2342.472 -2344.720 -2359.216
## 5 -2354.462 -2352.661 -2376.045 -2356.225 -2366.193 -2374.132 -2385.288
## 6 -2368.534 -2364.464 -2392.133 -2371.732 -2387.445 -2387.218 -2398.974
## 7 -2381.415 -2375.217 -2398.950 -2392.963 -2384.183 -2413.721 -2426.488
## 8 -2401.464        NA        NA -2385.815 -2404.950 -2434.531 -2434.990
## 9 -2404.102        NA        NA -2418.305 -2428.351 -2431.113 -2447.286
## 
## Top 3 models based on the BIC criterion:
##     EEE,3     EEE,4     VVE,2 
## -2314.386 -2320.207 -2320.433